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Abstract. We present the formula for angular distribution of 
integral flux of conventional (ir, K) muons deep under water 
taking into account the sphericity of the atmosphere and fluc- 
tuations of muon energy losses. The accuracy of this formula 
for various sea level muon spectra is discussed. The possi- 
bility of reconstructing two parameters of sea level spectrum 
by fitting measured underwater angular intensity is shown for 
Baikal Neutrino Telescope NT-36 experimental data. 



1 Introduction 

The knowledge of expected angular distribution of integral 
flux of atmospheric muons deep underwater is of interest not 
only for cosmic ray physics but also for the estimation of the 
possible background for neutrino detection and at last for a 
test of the correctness of underwater telescope data interpre- 
tation using the natural flux of atmospheric muons as cali- 
bration source. The last item frequently implies the estima- 
tion with an appropriate accuracy (e.g., better than 5 % for a 
given sea level spectrum) the underwater integral muon flux 
for various sets of depths, cutoff energies and angular bins 
especially for telescopes of big spacial dimensions. 

Up to now the presentation of the results of calculations 
of muon propagation through thick layers of water both for 
parent muon sea level spectra (especially for angular depen- 
dence taking into account the sphericity of atmosphere) and 
for underwater angular flux has not been quite convenient 
when applied to concrete underwater arrays. In addition, a 
part of numerical results is available only in data tables (of- 
ten insufficient for accurate interpolation) and figures. The 
possibility of direct implementation of Monte Carlo methods 
depends on the availability of corresponding codes and usu- 
ally assumes rather long computations and accurate choice 
of the grid for simulation parameters to avoid big systematic 
errors. Therefore, there remains the necessity of analytical 
expressions for underwater muon integral flux. In addition, 



the possibility of reconstructing the parameters of a sea level 
spectrum by fitting measured underwater flux in the case of 
their direct relation looks rather attractive. 

In this paper we present rather simple method allowing one 
to analytically calculate the angular distribution of integral 
muon flux deep under water for cutoff energies (1-10 4 ) GeV 
and slant depths of (1-16) km for conventional (ir, K) sea 
level atmospheric muon spectra fitted by means of five pa- 
rameters. The fluctuations of muon energy losses are taken 
into account. 

The possibility of reconstructing two parameters of sea 
level spectrum by fitting measured underwater angular inten- 
sity is shown for Baikal Neutrino Telescope NT-36 experi- 
mental data. 



2 Basic formulas 

According to the approach of work (Klimushin et al., 2001) 
the analytical expression for calculations of underwater an- 
gular integral flux above cutoff energy Ef for a slant depth 
R = h/ cos 8 seen at vertical depth h at zenith angle 8 and 
allowing for the fluctuations of energy loss is based on the 
relation 



F fl (>E f ,R,< 



F c i{>E f ,R, 8) 

c f {>E } ,R,ey 



(1) 
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where correction factor Ct is expressed, by definition, by the 
ratio of theoretical integral flux calculated in the continuous 
loss approximation to that calculated by exact Monte Carlo, 
and F c i(> Ef, R, 8) is the angular flux based on continuous 
energy losses. 

In principle, the correction factor Cj can be calculated us- 
ing known codes for muon propagation through matter. In 
this work we apply for this aim the MUM code described in 
work (Sokalski et al, 2001). 

The values of correction factors calculated for the same 
slant depth R at vertical direction and at zenith angle 8 differ 
weakly. It is illustrated in Fig. 1, where one can see that 
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Table 1. Coefficients Cij of the fitting formula (2) for correction factor calculated for vertical sea level spectrum given by expression (6) 
below. 
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C f (> E f ,R,0°) differs from C/(> E f , R, arccos/i/i?) 
maximum on 3.3 % for Ef >10 GeV at vertical depth h of 
1.15 km. It appears that with acceptable accuracy the correc- 
tion factor depends on slant depth R only, rather than on R 
and 9 separately. 

The dependencies of correction factor on Ef and R, cal- 
culated for sea level spectrum given by expression (6) below 
represent the set of rather smooth curves (shown in Fig. 1) 
and it is possible to approximate this factor by formula 



C f (> Ef, R, 9) = c v lo Sio E f) Rl 



(2) 



=0 j=0 



Here cut-off energy Ef is expressed in (GeV) and slant depth 
R is in (km) with the coefficients dj collected in Table 1. 
When using (2) for cutoff energies Ef <10 GeV one should 
substitute value of Ef =10 GeV. 

Formula (2) can be applied for any geometrical shape of 
the surface. Right hand side of (2) depends on 9 because, 
generally, R = R{9). So, in the particular case of a flat sur- 
face the angular dependence of the correction factor appears, 
in our approximation, only through the relation R — hJ cos 9 
(where h is a vertical depth). 

The accuracy of formula (2) for £/=(l-100) GeV is bet- 
ter than ±2 % for slant depths R as large as 22 km and is 
not worse than ±3 % for Ef=\ TeV up to i?=17 km and for 
Ef =10 TeV up to i?=15 km. Fig. 1 shows that for Ef < 100 
GeV the total energy loss may be treated as quasi-continuous 
(at level of C/ > 0.9) only for slant depths R < 2.5 km but 
for Ef =10 TeV the fluctuations should be taken into account 
at level of 15 % already for slant depth as small as R=l km. 

The dependence of correcton factor C/ on different sea- 
level vertical spectra is illustrated by Fig. 2. The correction 
factors calculated for £7/ =10 GeV using sea level spectrum 
(6) with spectral index 7 of 2.5 and 3.0 (instead of 2.72) differ 
more than on a factor of 2 starting from slant depth of i?=12 
km. Nevertheless, the values of C/ calculated using sea level 
spectra having 7=(2.65-2.78) are already within ±5 % corri- 
dor. For Ef=\ TeV this corridor is larger on 2 %. This fact 
results in the possibility to extrapolate the parametrization (2) 
based on sea level spectrum having 7=2.72 to other spectra 
at least up to slant depths of (12-13) km without introduction 
of additional spectral corrections. 

The angular flux F c i (> Ef,R, 9) based on effective linear 
continuous energy losses a + (3E having 2 slopes, is calcu- 
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Fig. 1. Correction factor Cf as a function of slant depth R in pure 
water. The results obtained using sea level spectrum defined by 
expression (6) are given. Solid curves correspond to numerical cal- 
culations for vertical case 6 = 0°. Dashed curves describe the cor- 
rection factor computed at vertical depth h of 1.15 km for various 
zenith angles as a function of slant depth defined by R = h/ cos 0. 
Both solid and dashed curves are shown for four values of cut-off 
energy Ef. 10 GeV, 100 GeV, 1 TeV, and 10 TeV, from top to bot- 
tom. 



lated by the following rule: 

F cl (>E f ,R,9) = 

Fd{>Ef,R,6\a u fo) for R<R 12 , 

F c i{>E l2 ,{R-R l2 ),9- a 2 ,fo) for R>R 12 . 

Here E\ 2 is the energy in the point of slope change from 
(ai, (3\) to (a 2 ,(3 2 ) and R 12 is the muon path from the en- 
ergy £12 till Ef which is given by 



(3) 



R12 = ^- In 

Pi 



«i + Ei 2 (ix 
ai + Effa. 



The formula for integral muon angular flux in the assump- 
tion of linear continuous energy losses is as follows: 

F cl (>E f ,R,6; a,0) = — — 
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Fig. 2. Correction factor C; as a function of spectral index 7 of 
sea level spectrum for various depths in water for vertical direction. 
The distributions for cut-off energy Ef =10 GeV are given. Solid 
curves correspond to numerical computations using sea level spec- 
trum defined by Eq. (6) with varying spectral index 7. Open circles 
correspond to numerical computations using VZK sea level spec- 
trum (Volkova et al, 1979), closed circles - Gaisser's sea level spec- 
trum (Gaisser, 1990), squares - MACRO (Ambrosio et al, 1995) 
sea level spectrum. All distibutions are shown for the following 
values of vertical depth in pure water: 1.15 km (a), 3 km (b), 5 km 
(c), 7 km (d), 9 km (e), 11 km (f), 13 km (g), 15 km (h), 17 km (i), 
and 21 km (j), from top to bottom. 
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where subscript i stands over both pion (tt) and kaon (K) 
terms and 



yi = ^(l-e-' 3R )+E c o r(6)e-e R , 

00 / n \ —1 

S(,, 7 ) = l + £n!z" I](7 + .7) = 1 + 
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When using expression (4) for slant depths R > i?i 2 one 
must substitute R — > (R — R12) and Ef — > £i 2 and use the 
values (a 2 , /3 2 ) for a loss description. For slant depths R < 
R12 the use of (4) remains unchangeable and the loss values 
are expressed by (ai , (3\). This algorithm may be extended to 
computations with any number of slopes of the energy losses. 

The 5 parameters (D ^ D 0k , E^(0°), E^ r K (0°),j) are 
those of the differential sea level muon spectrum, for which 



we use the following parametrization: 

-7 D °i 



D(E ,8) = E I 



° J^ K 1 + E o/ E S r A V 



(5) 



where 7 is a spectral index and E^ (9) have approximate 
sense of critical energies of pions and kaons for given zenith 
angle and E^ K {0°) are those for vertical direction. The 
corresponding angular distrubution should be introduced us- 
ing an analytical description of effective cosine cos 6* taking 
into account the sphericity of atmosphere. It should be noted 
that the description of underwater angular flux with the 5 pa- 
rameters of a sea level spectrum gives the possibility of their 
direct best fit using the experimental underwater distribution. 

Flux value in (4) is expressed in units of (cm _2 s _1 sr _1 ) 
and all energies are in (GeV), slant depth R in units of 
(gem -2 ), loss terms a and j3 in units of (10~ 3 GeVcm 2 g _1 ) 
and (10 _6 cm 2 g _1 ), correspondingly. For the description of 
effective linear continuous energy losses we use the follow- 
ing values of parameters when substituting in (3): (ai=2.67, 
/3i=3.40) and (a 2 =-6.5, /3 2 =3.66) with £i 2 =35.3 TeV. 

To examine the angular behaviour of a flux given by the 
formula (1) by means of the comparison with numerical cal- 
culations we used the following parameters of the sea level 
muon spectrum: 

D 07r = 0.175, D 0k = 6.475 x 1(T 3 , 
Efc (0°) = 103 GeV, E^ r K (0°) = 810 GeV, 7 = 2.72 . 

These values have been chosen according to splines com- 
puted in this work via the data tables kindly given us by au- 
thors of Ref. (Misaki et al., 1999). When checking the values 
of fit spectrum for cos #=(0.05-1 .0) we realized that the stan- 
dard description of effective cosine (with geometry of spher- 
ical atmosphere and with definite value of effective height of 
muon generation) is not enough and one should introduce an 
additional correction S(8) leading to (10-20)% increase of 
effective cosine value for cos 9 < 0. 1 . The reason of an ap- 
pearing of this correction is that the concept of an effective 
generation height is approximate one. It fails at large zenith 
angles where the real geometrical size of the generation re- 
gion becomes very large. 

The resulting fit of angular sea level spectrum in units 
of (cm _2 s _1 sr _1 Gey _1 ) is given by 

D(E ,6) = 0.175£(7 2 ' 72 

/ 

1 n:-!7 i 

(6) 
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with modified effective cosine expressed by 

cos0" = S(6) cos 9*, 



(7) 



where cos 9* is derived from spherical atmosphere geometry 
and is given by the polynomial fit: 



cos 9* = ^ a cos 1 9, 



(8) 



t=0 
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Fig. 3. Effective cosine as a function of zenith angle. Curve (a) 
is geometrical effective cosine cos 9* given by Eq. (8). Curve 
(b) is effective cosine cos#" with the correction and is given by 
Eq. (7). Curves (c) and (d) represent the ratio cos 9/ cos 6* and 
cos 9/ cos 9**, correspondingly. 

with the coefficients of the decomposition assembled in Ta- 
ble 2. The accuracy of (8) is much better than 0.3 % except 
the region cos #=(0.3-0.38) where it may reach the value of 
0.7 %. Note that for cos 9 > 0.4 the influence of the curva- 
ture of real atmosphere is less than 4 % but for cos 9 < 0. 1 it 
is greater than 40 % (Fig. 3). 

S{9) is the correction which is given for sec 9 < 20 by 

S{9) = 0.986 + 0.014 sec 9. (9) 

Correspondingly, for critical energies in expression (6) one 
should use cos 9** instead of cos 9*. 

The energy region, inside which the deviation of angular 
spectrum given by Eq. (6) from parent one is less than 5 %, 
is shifted from (0.3-200) TeV for cos 0=1.0 to (1.5-300) TeV 
for cos 0=0.05. The sea level spectrum given by (6) is valid 
only below the knee (Eq ~ 300 TeV) of primary cosmic ray 
spectrum. 

3 Comparison with numerical calculations 

The examination of (4) showed rather quick convergence of 
series S(z,j) with increase of R and Ef. Therefore, for 
the accuracy of F cl computation better than 0.1 % it is quite 
enough to take only four first terms of this series (up to z 3 ) 
for all values R > 1 km and Ef in (1-10 4 ) GeV. Even using 
the two terms leads to the accuracy of 1.3 % for (i?=1.15 km, 
E f =l GeV) and <0.5 % for (R > 2.5 km, E f > 1 GeV). 

Fig. 4 shows the comparison of underwater angular inte- 
gral fluxes allowing for loss fluctuations at different basic 
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Fig. 4. Underwater integral muon flux allowing for loss fluctuations 
as a function of zenith angle at different vertical depths. Four pic- 
tures are shown for various cut-off energies Ej: 10 GeV (a), 100 
GeV (b), 1 TeV (c), and 10 TeV (d), correspondingly. Four curves 
at each picture correspond to vertical depths h: 1.15 km, 1.61 km, 
2.0 km, and 3.0 km, from top to bottom. Solid curves result from 
numerical computations using the sea level spectrum based on data 
tables from (Misaki et ah, 1999) and MUM code of muon propaga- 
tion. Dotted curves result from analytical expression (1) using the 
sea level spectrum (6). 

depths h (of location of existing and planned telescopes) cal- 
culated both numerically using MUM code (Sokalski et al., 
2001) for parent sea level spectrum and analytically (1) for 
the spectrum given by (6). 

We realized that the error given by formula (1) for all men- 
tioned sea level spectra is within the corridor of ±(4-6) % for 
all cutoff energies £/=(l-10 3 ) GeV and slant depths R=(l- 
16) km (corresponding angle is expressed by cos0 = h/R 
for a given vertical depth h). This is proved for h in a range 
(1-3) km. For bigger cutoffs of £7/=(l-10) TeV the corridor 
of errors is ±(5-7) % for i?=(l-13)km. Note that for the 
sea level spectrum (6), just used for Cf parametrization, the 
errors are smaller on 2 %. 

The accuracy of the parametrization, used for the correc- 
tion factor as a function of Ef and slant depth R is rather 
high and is about ±5 % for all angles and kinds of the sea 
level spectrum (assuming that the spectral index 7 is approx- 
imately within (2.65-2.78)) (Fig. 5). It results in the possi- 
bility to use it for an estimating numerically from various sea 
level spectra the value of an angular integral flux allowing 
for fluctuations of losses without direct Monte Carlo simula- 
tions. 

Note that the expression (1) may be directly used for an ice 
after the substitution R — > R/p, with p being the ice density, 
and, with an additional error of ~ 2 %, for sea water. In spite 
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Table 2. Coefficients c, of the fitting formula (8) for effective cosine with the maximum relative errors. 



cose 


CO 


Cl 


C2 


C3 


C4 


Max.err,% 


CH-0.002 


0.11137 














0.004 


0.002=0.2 


0.11148 


-0.03427 


5.2053 


-14.197 


16.138 


0.3 


0.2=0.8 


0.06714 


0.71578 


0.42377 


-0.19634 


-0.021145 


0.7 



described in Sec. 1 we have realized that: 
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Fig. 5. The accuracy of formula (1) as function of cutoff energy Ef 
for different vertical sea level spectra. Four pictures are shown for 
various vertical depths h: 2 km (a), 5 km (b), 10 km (c), and 15 
km (d), correspondingly. Thick solid curves correspond to compar- 
ison with numerical computations using sea level muon spectrum 
based on data tables from (Misaki et ah, 1999) and MUM code of 
muon propagation. Thin solid curves - sea level spectrum defined 
by Eq. (6), dashed - VZK sea level spectrum (Volkova et ai, 1979), 
dash-dotted - Gaisser's sea level spectrum (Gaisser, 1990), dotted - 
LVD (Aglietta et al, 1999) sea level spectrum. 



of seeming complexity of the formulas (1), (3) and (4) they 
may be easily programmed. 

The validity of this analytical expression with an accuracy 
of ±(5-7) % for £/=(10 3 -10 4 ) GeV and slant depths of (1- 
12) km gives also the possibility of estimation the angular 
underwater differential spectrum (by means of numerical dif- 
ferentiation) with error smaller than ±(6-8) % for energies of 
(30-5 xlO 3 ) GeV. 



4 Parametrization of atmospheric muon angular flux 
using underwater data 

When reconstructing the parameters of sea level spectrum 
defined by Eq. (5) by fitting with MINUIT least square 
method the corresponding underwater angular intensity ex- 
pected at vertical depth h= 1 . 1 5 km and expressed by formula 



(i) it is possible to reconstruct two parameters (-Do, > l) of 
sea level spectrum when angular bins corresponding to 
slant depth R >6 km are involved 

(ii) the reconstruction of third parameter Dq k is formally 
possible only using angular bins corresponding to slant 
depths R >15km where neutrino induced intensity 
should be taken into account. 

For checking the same procedure using experimental re- 
sults we have examined the data sample with NT-36 (1993) 
unfolded experimental angular intensity published by Baikal 
Collaboration in Ref. (Belolaptikov et ai, 1997) for vertical 
depth of ft=1.15 km. The cutoff energy value was taken as 
.E/=10GeV The whole data sample corresponds to 44 an- 
gular bins A cos 6=0.02 (cos 6>=(0. 13-0.99)) with maximum 
slant depth i?=8.8km. The mean muon energy at the sea 
level corresponding to this angular range is E=(0.6-15) TeV 
Only statistical errors have been taken into account. The fol- 
lowing 3 parameters of sea level spectrum were taken accord- 
ing to expressions (5) and (6): 

D 0k = 0.037A>., 
Ef£ (0°) = 103 GeV, E^ r K (0°) = 810 GeV. 

The results of reconstructing of two free parameters (D 0jr , 7) 
of sea level spectrum are as follows. 

(i) For a range of zenith angles within cos 8=(0. 17-0.99) 
we have obtained formally (£> 07r = 0.26, 7 = 2.79). It 
is illustrated by Fig. 6. In spite of this result coincides 
with MACRO (Ambrosio et ai, 1995) and LVD (Agli- 
etta et ai, 1999) best fits, its confidence level (CL) is 
close to 0. The artificial increase of errors in 3 times 
due to additional systematic errors leads to (Dq^ = 
0.17,7 = 2.73) with CL=87%. 

(ii) For vertical directions with cos 8=(0. 6 1-0.99) the re- 
constructed sea level spectrum is extremely steep with 
(Do, =1.0,7 = 3.0) and CL=0.5 % but the increase of 
errors in 2 times results in (Z? „ = 0.19, 7 = 2.74) with 
CL=40 %. 

(iii) For horizontal directions with cos 8=(0. 13-0.61) the 
reconstructed sea level spectrum is flat, as (D 0jr = 
0.1,7 = 2.65) with CL=70%, and for cos(9=(0.17- 
0.61) as (D 07i = 0.12,7 = 2.68) with CL=40%. The 
result of this best fit is shown in Fig. 7. 
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Fig. 6. Zenith angle distribution of the muon intensity at vertical 
depth of 1.15 km. Experimental points - NT-36 data (Belolaptikov 
et al, 1997). Solid curve results from the best fit of 2 parameters 
(Do T = 0.26, 7 = 2.79) of sea level spectrum. Dashed curve re- 
sults from analytical expression (1) using the sea level spectrum (6) 
and is consistent with experimental data also with CL=0. Only sta- 
tistical errors have been taken into account. 



Fig. 7. Zenith angle distribution of the muon intensity at vertical 
depth of 1.15 km for horizontal directions. Experimental points - 
NT-36 data (Belolaptikov et al, 1997). Solid curve results from 
the best fit of 2 parameters (-Do x = 0.12, 7 = 2.68) of sea level 
spectrum. Dashed curve results from the best fit of 2 parameters 
(Do, = 0.26, 7 = 2.79) of sea level spectrum using a whole an- 
gular range. Only statistical errors have been taken into account. 



It should be pointed out that the implementation of Gaisser's 
set of 3 parameters (D 0k , (0°), E^ K (0°)) (Gaisser, 1990) 
gives almost the same results of reconstructing of (Dq^ , 7), 
as well as when using the recalculated depth-intensity curve. 
The fact that sea level spectrum changes the slope from ver- 
tical directions to horizontal ones may be explained probably 
by unproper taking into account the muon bundles when un- 
folding the measured intensity. 

5 Conclusions 

The analytical expression presented in this work allows to es- 
timate for fluctuating losses the integral flux of atmospheric 
muons in pure water expected for different zenith angles, 
cos #=(0.05-1.0), at various vertical depths at least of h=(l- 
3) km for different parametrizations of the sea level muon 
spectra. The errors of this expression are estimated to be 
smaller than ±(4-6) % for cutoff energies ranged in Ef=(l- 
10 3 ) GeV and slant depths in h/ cos #=(1-16) km. The main 
advantage of the presented formula consists in the possibil- 
ity of the direct best fit of at least 2 parameters of parent sea 
level spectrum using angular distribution of underwater inte- 
gral flux measured experimentally at a given vertical depth. 

The fitted sea level spectrum for NT-36 data is too steep 
for vertical directions (7=3.0) and flat for horizontal ones 
(7=2.65-2.68). It leads to the necessity of proper introduc- 
ing of systematic errors mainly resulted from muon bundles. 



The artificial increase of statistical errors in 2-3 times results 
in sea level spectra closer to (Klimushin et al, 2001) and 
(Gaisser, 1990). 

The proposed method may be adapted to estimations in 
rock after corresponding description of the correction factor 
and continuous effective losses. 
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